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Abstract 

This paper describes and analyzes the spatial spread of tuberculosis (TB) on complex metapopulation, 
that is, networks of populations connected by migratory flows whose configurations are described in terms 
of connectivity distribution of nodes (patches) and the conditional probabilities of connections among 
classes of nodes sharing the same degree. The migration and transmission processes occur simultaneously. 
For uncorrelated networks under the assumption of standard incidence transmission, we compute the 
disease-free equilibrium and the basic reproduction number, and show that the disease-free equilibrium 
is locally asymptotically stable. Moreover, for uncorrelated networks and under assumption of simple 
mass action transmission, we give a necessary and sufficient conditions for the instability of the disease- 
free equilibrium. The existence of endemic equilibria is also discussed. Finally, the prevalence of the TB 
infection across the metapopulation as a function of the path connectivity is studied using numerical 
simulations. 
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AMS Classification: 34A34, 34D23, 34D40, 92D30 



1 



1 Introduction 



Despite significant advances in medical science, infectious diseases continue to impact human populations 
in many parts of the world. Tuberculosis (abbreviated as TB for tubercle bacillus) is a common deadly 
infectious disease caused mainly by Mycobacterium tuberculosis. It basically attacks the lungs (pulmonary 
TB), but can also affect the central nervous system, circulatory system, the genital- urinary system, bones, 
joints and even the skin. Tuberculosis can spread through cough, sneeze, speak, kiss or spit from active 
pulmonary TB persons. It can also spread through use of an infected person's unsterilized eating utensils 
and in rare cases a pregnant woman with active TB can infect her foetus (vertical transmission) [1,2]. 
Transmission can only occur from people with active TB but not latent TB. This transmission from one 
person to another depends upon the number of infectious droplets expelled by a carrier, the effectiveness 
of ventilation, duration of the exposure and virulence of the MTB strain. The chain of transmission can 
therefore be broken by isolating patients with active disease and starting effective anti-tuberculosis therapy 
[1-5]. At present, about 95% of the estimated 8 million new cases of TB occurring each year are in developing 
countries, where 80% occur among people between the ages of 15-59 years [1]. In sub-Saharan Africa, TB is 
the leading cause of mortality and in developing countries, it accounts for an estimated 2 million deaths which 
accounts for a quarter of avoidable adult deaths [1]. It is known that factors such as endogenous reactivation, 
emergence of multi-drug resistant TB, and increase in HIV incidence in the recent years call for improved 
control strategies for TB. A full understanding of the effectiveness of treatment and control strategics within 
different regions of the world is still needed. It is worth emphasizing that mathematical analysis of biomedical 
and disease transmission models can contribute to the understanding of the mechanisms of those processes 
and to design potential therapies (see [6-9] and references therein). A number of theoretical studies have 
been carried out on the mathematical modeling of TB transmission dynamics [3-9,38,39]. 

However, the analysis of the spread of infectious diseases on complex networks has become a central 
issue in modern epidemiology [10] and, indeed, it was one of the main motivations for the development of 
percolation theory [11]. While the initial approach was focusscd on local contact networks [12-16], that is, 
social networks within single populations (cities, urban areas), a new approach has been recently introduced 
for dealing with the spread of diseases in ensembles of (local) populations with a complex spatial arrangement 
and connected by the migrations. Such sets of connected populations living in a patchy environment are 
called metapopulations in ecology, and their study began in 1967 with the theory of island biogeography 
[17]. 

Unfortunately, when considering dispersal models, there is an approach based on the metapopulation 
concept. The population is subdivided into a number of discrete patches which are supposed to be well 
mixed. Then, in each patch the population is subdivided into compartments corresponding to different 
epidemic status. This leads to a multi-patch, multi-compartment system. At this point two formulations 
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are possible. 

The first one assumes that an infective in one patch can infect susceptible individuals in another patch. 
This assumption gives rise to a family of models which have been well studied [18,19]. This formulation 
assumes that there is a spatial coupling between patches, but that individuals (vectors or hosts) do not 
migrate between patches. They make short 'visits' from their home patches to other patches. 

The second one considers migration of individuals between patches. The infection does not take place 
during the migration process. The situation is that of a directed graph, where the vertices represent the 
patches and the arcs represent the links between patches. Recently, there has been increased interest in these 
deterministic metapopulation disease models. For instance, in some recent models of epidemic spreading, the 
location of the patches in space is treated explicitly (without taking into account the number of connections 
k (degree) that any given patch in the network may have) thanks to the increasing of computational power 
(see for instance [20, 21]). In Refs. [16, 22, 23], however, an alternative approach based on the formalism 
used in the statistical mechanics of complex networks is presented. Under this approach, the structure of 
the spatial network of patches (nodes) is encapsulated by means of the connectivity (degree) distribution 
p{k) defined as the probability that a randomly chosen patch has connectivity k. In contrast, in [24, 25], 
the authors consider reaction diffusion processes to take place simultaneously, which turns out to be correct 
assumption for a suitable continuous-time formulation of metapopulation models for the spread of infectious 
diseases. 

In this paper, we consider the spread of TB on complex metapopulations, that is, networks of populations 
connected by migratory flows whose configurations are described in terms of the conditional probabilities of 
connections among classes of nodes sharing the same degree. For uncorrelated networks under the assumption 
of standard incidence [37] (or frequency-dependent) transmission, we compute the disease-free equilibrium 
and the basic reproduction number and show that the disease-free equilibrium is locally asymptotically 
stable. Moreover, for uncorrelated networks and under assumption of simple mass action [37] (or density- 
dependent) transmission, we give a necessary and sufficient conditions for the instability of the disease-free 
equilibrium. We find that there exists a more precise bound of the largest eigenvalue of the Jacobian matrix 
of the system around the disease- free equilibrium. This condition says that, for fixed values of the migration 
rates of latently-infected and infectious individuals, a high enough density of individuals and/or large enough 
maximum connectivity in the metapopulation guarantee the instability of the disease-free equilibrium and, 
hence, TB spread. In the limit of infinite networks with bounded average degree, this condition implies 
the existence of a TB threshold for any distribution with large value. The existence of endemic equilibria 
is also discussed. Additionally, through numerical simulations, the forecasted prevalence of the infection is 
not constant but increases with the patch connectivity. Interestingly, close the epidemic threshold, there 
are always patches with low connectivities where TB is not able to progress unless infectious individuals 
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arrive from (crowded) patches with higher connectivities. Comparing to existing resuhs in the literature, 
our work treats a specific disease which is not the case in Refs. [24, 25, 26]. We point out that in Refs. 
[24, 25, 26], the authors have neglected some important epidemiological features of the propagation of a 
disease such as births, natural mortality, mortality due to the disease, natural recovery and the basic models 
studied arc of dimension 2 which arc very simple. In addition, the authors have supposed that the total 
population is constant which is not always the case. Our basic model is of dimension 4 and incorporates the 
essential biological and epidemiological features of TB such as birth, mortality due to the disease, slow and 
fast progression, effective chemoprophylaxis of latently-infected individuals, natural recovery and treatment 
of infectious, relapse from the disease and re-infection after recovery. Also in our model, the total population 
is not constant. It is our view fact that this study represents the first work that provides an in-depth the 
spread of TB on complex metapopulation using a degree of distribution and conditional probabilities. 

2 A TB metapopulation model 
2.1 The model 

Wc consider the spread of TB in heterogeneous metapopulations. The model consists of ?? patches rep- 
resenting n different degree of connectivities. We assume that the architecture of the network of patches 
(nodes) where local populations live is mathematically encoded by means of the connectivity (degree) dis- 
tribution p(fc), defined as the probability that a randomly chosen patch has degree k. At any given time, in 
each patch, an individual is in one of the following states: susceptible, latently infected (exposed to TB but 
not infectious), infectious (has active TB) and recovered. These states arc average number (density) of ps^k, 
PE,k , pi,k and pR^k in the patches of connectivity k, respectively. The total variable population size at time 
t is given by, 

Pk{t) = ps,k{t) + PEM{t) + piMt) + PbA^)- (1) 

It is assumed that births are recruited into the population at per capita rate A. The transmission of 

Mycobacterium tuberculosis occurs following adequate contacts between a susceptible and infectious in each 

sub-population. The rate at which susceptible are infected is i^ELhR—. for standard incidence (or frequency- 

Pk 

dependent) transmission and (3 pi.kPs.k for simple mass action (or density-dependent) transmission, where (3 
is the effective contact rate of infectious that is sufficient to transmit infection to susceptible (it also denotes 
how contagious of the disease is). On adequate contacts with active individuals, a susceptible individual 
becomes infected but not yet infectious. A fraction q of newly infected individuals is assumed to undergo a 
fast progression directly to the infectious class, while the remainder is latently infected and enter the latent 
class. Latently infected individuals are assumed to acquire some immunity as a result of the infection, which 
reduces the risk of subsequent infection but does not fully prevent it. We assume that chemoprophylaxis 
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of latently infected individuals reduces their reactivation at a rate 9 and that the initiation of therapeutics 

immediately remove individuals from active status and place them into a latent state. This last assumption 

is realistic. Indeed, the classic works of Jindani et al. [40] showed that a bactericidal treatment reduced 

the number of bacilli 20 times during the first two days and about 200 times during the 12 days. After two 

weeks of treatment, the sputum of a patient contain on average 1000 times less bacilli before treatment, a 

number generally too low to be detected on direct examination. Latently infected individuals who received 

successful chemoprophylaxis can recover at a constant rate t] and enter the recovered class R. Those who did 

not received effective chemoprophylaxis progress to active TB at a rate q;(1 — 9) where a is the rate at which 

latently infected individuals become infectious (this value is connected with the average time of incubations) . 

Once in active stage of the disease, due to their own immunity, an individual may recover naturally and 

will move in the class of latently infected at rate 7. Also, after a therapy of treatment, infectious can 

spontaneously recover from the disease and will enter the recovery class R at rate 5. As suggested by Styblo 

[41], recovered individuals can only have a partial immunity. Hence, they can relapse from the disease with 

a constant rate ^ and enter the infectious class /, in the same time some of them can be re-infected at a rate 

(1 — ^)/3^^ and enter class E. The rate for non-disease related death is /i, thus, 1//^ is the average lifetime. 
Pk 

Infectious have addition death rate due to disease with a rate d. 

This model description is summarized in the flow diagram given below. 

According to the derivation in [24, 25] of the continuous-time formulation for the progress of diseases on 
metapopulations and assuming non-limited or frequency-dependent transmission, the equations governing 
the dynamics of TB propagation are 

PE,k = /?(l-g)^^^+/3(l-0^^^^+7P/,.-[M + 'y + «(l-eW 
Pk Pk 

- DEPE,k + kDE E P{k'\k)^, 

k' 1^ 



(2) 



Pi,k = Pq^^^^^^^^ + a{l - 9)pE.k - {P' + d + ^ + 5)pi,k + i PF!^,k 
Pk 

~ Dipi,k + kDiY.P{k'\k)^, 

k' ^ 

PR,k = (1 - ^jP^^hfM. + rjpE,k + Spi,k -ip + O PB..k - Dn pn,k + kDn E P{k'\k)^'''''' 



Pk ' ' ' ' k' k' 

where k is the degree of the patches where local population live (k = ki, . . . , kmax), and P(k'\k) is the 
conditional probability that a patch of degree k has a connection to a patch of degree k'. As in classical 
reaction-diffusion processes, Eq. (2) expresses the time variation of susceptible, latently infected individuals, 
recovered individuals and infectious as the sum of two independent contributions: reaction and diffusion. In 
particular, the diffusion term includes the outflow of individuals (diffusing particles) from patches of degree 



k and the inflow of migratory individuals from the nearest patches of degree k' . For the sake of brevity, in 
the sequel we consider strictly positive diffusion rates {Dg, De, Dj, Dj^ > 0). 

For limited or frequency-dependent transmission model, we simple replace in Eq. (2) the transmission 

, n Pl,kPS,k , a 

term P by ppi^kPs,k- 

Pk 

2.2 Positively-invariant set 

Notice that, since births and deaths are considered in model (2), the total number of individuals is not 
constant at the mctapopulation level. More precisely, multiplying equations in system (2) by p(fc), and 
summing over all k, we have the following differential equations for p5, pi and p^f, the average number 
of susceptible, latently infected, infectious and recovered individuals per path at time t, respectively. 



PS 



PE 



A _ /?^p(fc)^£:^ - ^,p, _ DsPs + E E kp{k)P{k'\kf-fL, 

P{l-q)Y. p{kf-i^^ + /3(1 - E ^^^^ +^pi-[p + r^ + a{l- 6)]pE 
k k 

- DePe + DEY,Y.^p{k)P{k'\k) ^ , 

A; k' 

PI = PqY^p{k)^i^^+ail-e)pE-ip + d + ^ + 5)pi (3) 



-Dipi +DjY,Y1 kp{k)P{k'\k)P^ +ipR, 

k k' 



PR = -/^(i - oy]p(^) ^^''^^^''' +VPE + Spi - (h + Opb^ 

k p^ 

- Dnpn + ^i? E E kp{k)P{k'\k)^, 

k k' 

where Pj{t) = ^^p(fc)pj_fc, j = S,E,I,R. Now, since the number of hnks emanating from nodes of degree 
fc 

k to nodes of degree k' must be equal to the number of links emanating from nodes of degree k' to nodes of 
degree k in non-directed graphs, we have the following relationship between p{k) and P{k'\k) [27]: 

kP{k'\k)p{k) = k'P{k\k')p{k'). (4) 

Using this restriction and the fact that 'Y^P{k\k') ~ 1 after changing the order of summations, Eq. (3) 

k 
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becomes 



PS = J^-P2^P[k) MPs, 



Pk 

(5) 



k Pk ^ Pk 



Pi = (3qY,P{k) ^^'^^^'^ +a{l - 9)pE - {p + d + + S}pi + ^ pr, 
k Pk 



PR 



= - E P^^) ^^"^^ + PE + SpI-{^, + OPR■ 



k 



Adding the expressions in the right-hand side of the equations in system (5) yields 

^^A-pp~dpi. (6) 

From the above equation, one can deduce that -j- < A — fip. Thus, ^ < if p > — . Since < A — pp, it 

at at n at 

can be shown that using a standard comparison theorem [28] , that 

pW<P(0)e-^* + -(l-e-^*). 
M 

If p(0) < -, then p{t) < -. 

P P 

Hence, all feasible solutions of components of system (5) enters the region: 

(^{ps,pe,pi,pb) eRio^ pit)<j^Y (7) 

Thus, it follows from Eq. (7) that all possible solutions of system (5) will enter the region il. Hence, the 
region fi, of biological interest, is positively-invariant under the flow induced by system (5). Further, it can 
be shown using the theory of permanence [26] that all solutions on the boundary of eventually enter the 
interior of fl. Furthermore, in fi, the usual existence, uniqueness and continuation results hold for system 
(5). Hence, system (5) is well posed mathematically and cpidemiologically and it is sufficient to consider the 
dynamics of the flow generated by system (5) in fl. The same conclusions on hold for the simple mass 
action (or density-dependent) model. 

For networks with a connectivity pattern defined by a set of conditional probabilities P(k'\k), we define 
the elements of the connectivity matrix C as 

Note that these elements are the average number of individuals that patches of degree k receive from 
neighboring patches of degree k' assuming that one individual leaves each of these patches by choosing at 
random one of the k' connections [9]. One should notice that, for those degrees k that are not present in 
the network, P(k'\k) — 0, Vfc'. Hereafter in the paper, when talking about degrees, we implicitly mean 
those degrees that are present in the network. Furthermore, the case with patches having all the same 
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connectivity is excluded from our considerations because, under the present approach, the model equations 
reduce to those of a single patch SEIR model. 

3 Uncor related networks 



In order to obtain analytical results about the TB metapopulation dynamics, we need to be precise about 
the form of P{k'\k). The easiest and usual assumption is to restrict ourselves to uncorrelated networks. In 
these networks, the degrees of the nodes at the ends of any given link are independent, that is, no degree- 
degree correlation between the connected nodes. In this case, we have that P{k'\k) ~ k'p{k')/{k) which 
corresponds to the degree distribution of nodes (patches) that arrive at by following a randomly chosen link 
[29]. 

3.1 Analysis of standard incidence (or frequency-dependent) model 

After replacing the expression of P{k'\k) into Eq.(2), one obtains the following equations for TB spread 
in metapopulations described by uncorrelated networks and limited transmission: 



Ps,k 



PE,k 



f. oPl.kPS.k p, / ^" ^ 

A - P MPs.fc - Ds Ps,k - -tjtPs , 

Pk \ [k) J 

^Pl,kPS,k , , .^PUkPRM , r I I n AM 

P(l - 9) + +P(1 - + IPi.k -[P- + V + a(l - 0)\pE.k 

Pk Pk 



Pl,k 



PR,k 



-De [ PE,k - J^PE 



p^PijkPsj^ + a(l - 9)pE.k - {p + d + ^ + 5)pi^k - Di \ pi.k - -nrPi 

Pk ' V (fc) . 

_ ^) Pl^kPR^k ^ ^^^^ ^ ^ _ + ^) J. _ f p _ JL 

Pk V («> 



(8) 



^PR,k, 



PR 



where (fc) = ^ kp{k) is the average network degree, 
fc 

In this form, it becomes clearer that the diffusion term is simply given by the difference between the 
outflow of susceptible, latently infected, infectious and recovered individuals in patches of connectivity fc, 
DsPs.k, DEPE,k , Dipi k, and Drpbm and the total inflow of susceptible, latently infected, infectious and 
recovered individuals across all their k connections, which is k times the average flow of individuals across 
a connection in the network, Dsps/{k), DEPE/{k) , Dipi/{k) and Dnpji/ {k). Note that this average flow 
across a connection does not depend on the degree k of the considered patch because we have assumed that 
the architecture of the metapopulation is described by an uncorrelated network. 
In these networks, the elements of the connectivity matrix C are simply 

kp{k') 



Ckk' — 



(k) 



(9) 



Clearly, C is a rank-one matrix and has the vector with components Vk — k as eigenvector of eigenvalue 
1. So, if there are n different degrees in the network, then the eigenvalues of this matrix arc A = 0, with 
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algebraic multiplicity n — I and A = 1 which is a simple eigenvalue. This fact will be used in the stability of 
equilibria of the model. To do this, we are going to 'vectorializc' system (8), using the following vectors of 



'5' = {Ps.ki, PS,k2i ■ ■ ■ 1 PS,kr,V , E = {pE.ki, PEM' • ■ • ' PE,k„)^ , I = (Pl.ki, PiM' • ' • ' Pl.k^Y , 

R = {PR,ki^PR,k2^---^PR.k„V, N = {pk,,Pk2, - ■ ■ ,Pk„V and I = (1, 1, . . . , 1)^. If X e M" is a vector, we 
denote by diag(X) the n x n matrix whose diagonal is given by the components of X. With these notations 
and conventions, system (8) becomes 

' S = AI - l3dmg{N)-^di&g{I)S - {p + Ds)S + DsCS, 

E = - g)diag(A^)-idiag(/)5 + /3(1 - S)'^i&g{N)~^Amg{I)R 

+ [/i + 7] + a{l -e) + De]E + DeCE, (10) 

/ = l3qAmg{N)-^A\s.g{I)S + a{l-e)E-{ji + d + -f + 5 + Di)I + DiCI + £,R, 

, R = -^l3{l~^)dia.g{N)-^diag{I)R + r]E + 6I-{p + ^ + DR)R + DRCR, 
where C is the connectivity matrix defined as in Eq. (9). 

We point out that in the case where the parameters /3, q, 7, p, S, 9, a, ^, 77 and d are not the same for all 
patches, they are replaced in system (10) by diagonal non-negative matrices and this does not change the 
fundamental structure of the system. 

3.1.1 Disease- free equilibrium (DFE) for generic networks 

The disease- free equilibrium of model system (2) are the solutions p'g p% and p'j ^ to the equations: 

A - pPlh^ _ ,,pO _ Dsp%^, + kDs Y: Pim'^ = 0, 
Pk fc, 



sPi.kPs,k , an ^^Pi.kPR^k ,0 r I 1 n AM 
^(1 - 1) 1) + - + lPi,k - [P + ri + a{l - e)]p% J, 



Pi 



Pi 



P%.k' 



-DEP%^k + kDE E P{k'\k)-^ = 0, 



(11) 



Pi 



Pl.kPR,k 



pI 



lPE,k + S pIu - (m + e) Pfl.fe - Dr + kDn J2 P{k'\k)^ = 



For the analysis of the infection's spread, the so-called disease-free equilibrium is particularly relevant. 
By definition, this is obtained by replacing pi^k = in Eq.(2), leading to an explicit expression for the 
number of susceptible individuals in patches with degree k that can be written as 

{p + Ds)p%,=A + DsJ2Ckk'Plk'- 



As P{k'\k) — 1, it follows that, for any generic network, one has 

k' 



Note that Eq. (6) at the disease-free equihbrium yields 



Then, the disease-free equilibrium is given by 



(12) 



P + Ds \ ' ^J. (k) ^ 

Equation (12) is also the disease free equilibrium for the simple mass action transmission model. 
3.1.2 Basic reproduction number and local stability of (DFE) 

The global behavior for this model crucially depends on the basic reproduction number, that is. the 
average number of secondary cases produced by a single infective individual which is introduced into an 
entirely susceptible population. System (8) has an evident equilibrium Qq — (S*", 0,0,0) with ~ p^g j. 
defined as in Eq. (12) and the zero vector of dimension n when there is no disease. We calculate the basic 
reproduction number, 7?.o, using the next generation approach, developed in Ref. [30]. 

Using the notations in Ref. [30], the matrices F and V , for the new infections and the remaining transfers, 
are, respectively, given by 



F = 



Fi 0" 





and V = 



AeIu-DeC -7^ 
-a (1 - e) /„ Ai /„ -DiC -iln 

-Vln -Sin AnJn - DrC 



where /„ is the identity matrix of dimension n, 



Fi = /3(1 - g)/„. Fa = /3gf/„, = [/i + ry + a (1 - 6*) + D^], Ai ^ p + d + ^ + 5 + Di, An^fi + ^ + Di 



Set 



where 



Also, let 



where 

V, = 



Fu = 



Fi 
F2 



F 



F12 = 



V 



Fii F12 
F21 F22 



F21 = [0, 0] and F22 = 0. 



Vi V2 



Ae In -DeC --fin 

-a{l-e)In AiIn-DiC 



V2 







, V3^[-T]In, -Sin] and Vi^ARln-DnC. 



We stress that since is a M-matrix and —V is stable , then one can deduce that V ^ > 0. 
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Now, we need to compute the inverse of the matrix V. To this end, suppose that the inverse matrix of 
V can be written in the fohowing form: 



Wii Wi2 

W21 W22 



where Wu and W22 are square matrices of dimension (2n x 2n) and {n x n), respectively. 
Observe that 



A B 




where A — Fn Wu and B ~ Fu W\2- Thus, the basic reproduction ratio is defined, fohowing [30], as the 
spectral radius of the next generation matrix, FV~"^\ 



Tlo = p{FV'^) = p{A) = p {Fn Wu) . 



(13) 



To compute the explicit expression of the basic reproduction number, we need to compute the inverse matrix 
of V. To this end, we need the following lemma stated above and proved in Appendix A. 



Lemma 1 ; Let N be a square block matrix of the following form: 

N 



N2 



where Ni and N4 are square matrices. 

If Ni and D^Ni- A^gA^f ^iVa are invertible, then the inverse matrix of N is given by 

'N^^ + N-^N2D-^N3N[^ -N-^N2D- 

Note that the matrix Vi has the form of the matrix A^ defined in Lemma 1 with A^i — Ae In — DeC, 
N2 = -7 In. N3 = -a{l - e) /„ and = Ai /„ - Di C. 

Note also that the matrix V has the form of the matrix A^ defined in Lemma 1 with A'l = Vi, A^2 = ^^2, 
A'a = ^3 and A4 = V4. So, if all the hypotheses in Lemma 1 arc satisfied for the matrices Vi and V4, then 
Lemma 1 can be used twice to compute V~^. 

Thus using Lemma 1, one can prove that V^'^ has the following form: 



v; 



Vii V12 

V21 V22 



where 



^11 = {AeIu-DeC)-' +-/ {AeIu -De Cy' V2U 
V12 = i{AEln-DEC)-^ 

V21 = a{l~0)V22iAEln~DEC)-\ 



22 



Ai In ~ DiC -aj(l-e) {Ae In - De Cy 
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From the above expressions, it appears that to compute the expHcit expressions of Vn, V12, V21 and V22, we 



need to compute the inverse matrices of 



Aj /„ ~DjC~a^{l~e) {Ae In - De Cy^] and {Ae In ~ De C). 



To do so, we shall used the following Lemma 2 stated below and proved in Appendix C. 

Lemma 2 ; Let G = U + XWZbeannxn invertible matrix. Suppose that the matrices U , W and 
+ ZU^^ X are invertible. Then, the inverse matrix of R is defined as 



(14) 



Using the above Lemma 2 and the fact that C" = C, Vm E N* , one can easily prove that 

' De ' 



{AeIu-De cy' = ^ 

Ae 



Ae-De 



C 



Ai /„ -DiC-oi^{\-B) {Ae /„ - De Cy^ 



-1 1 



a — b 



■C 



where 



Ai{fi + r] + DE) + a{l-e){fi + d + S + Dj) , , AeDj[^i + 'q + a{l - 9)] + ■y a{l - i 

a = : and — 



Af 



With this in mind, after some substitutions, one has: 



1^22 



V21 



12 



In + 



a — b 
ao In + bo C, 
a(l -61) 



C 



oAe 
oi /„ + bi C, 

1 



b\^ + a{)^_ffj\+aDE^ 
J^n H ; TTT Tz :tt^l^ 



{a - b)[^i + a{l ^ 9)] 



aAf 



^ b[^i + ail~9)]+aDE ^ 



{a ~ b)[^i + ail - 9)] 



0,2 In + h C, 

oAe + 7 a{l 
a Ae + 7 a{l 



■ ^a{l^9)AE[b[^i + a{l-9)]+aDE\ 

(a-5) [aAE+ia{l-9)] [^l + a[l - 9)Y 

'DE[^l + a(l - &)] [oAe+I a{l - 9)] 



[a Ae + 1 a{l - 9)] [^i + a{l-9y 



as In + h C, 
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where, 



ao 



bo 



ai 



hi 



a2 



1 ^ Ae 

a Ai{ji + 7^ + DE)+a{l-e){^i + d + 5 + Di)' 

b 



a {a --by 

a(l -6>) 
ciAe 

a(l - e) b[^JL + a{l-e)]+aDE 

qAe {a ^ b)[^i + a{l - e)] ' 

1 



as 
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aAs ' 

1 b[n + a{\~e)]+aDE 
qAe (a-5)[/i + a(l - 6I)] ' 

a^_E +7a(l - 6*) 



aAE + 7 a(l — 



aAl 



7 a (1 - 6*) As [6 [a* + q(1 - 61)] + uDe] + {a-b) De [a* + q(1 - 6I)] [a + 7 a(l - 9)] ' 



{a - b) [a Ae + -f a{l - e)][^i + a{l - 9)]^ 

This achieve the computation of V{~^. 

Now, we need to compute V~^. To this end, we need to prove the invertibihty of matrix D ~ V4 
V3 V-f^ ¥2- Simple substitutions show that: 

D = Vi-i7^^Vii+6^V2i), 

= [AR-Cir]a3 + 5ai)]Ir,-[DR + ^{rjb3 + Sbi)] C. 
Applying Lemma 2 one again, the inverse of D is given by 



1 



DR + i{iib:i + 5bi) 



where 



04 



[AR-S,{iia3 + 5ai)] 

CI4 In — 64 C. 

1 

[AR-£,{Tias. + 6ai)\' 
1 



[Ar ~i{vaz + 5 ai )] - [Dr + ^ (?? ^^3 + (5 61 )] 



C 



DR + Cinbs + Sbi) 



[AR-^{ija3 + 6ai)] [Ar - ^ [r] as + 6 ai)] - [Dr + ^ (r^bs + 6 bi)]' 

Since Vi and D arc invertible matrices, applying Lemma 1, after simple calculations we have: 

Wii = V,-' +V^-'V2D-^V3V^-\ 

>ii + e ^12^?"' {v Vn + S V21) V12 + ^ Vi2D~^ (7y Vu + S V22) 
¥21 + e V22D-^ (r; Vii + 6 V21) V22 + C V22D-^ (v V12 + S T^22). 
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At this stage, we need to compute the expression of A. Note that A can be written as follows: 
A = FnWn, 

Fi [V21 + CV22D-^ {l^Vn + SV21)] Fi [V22+iV22D-^{T^Vi2+5V22)] 

F2 [V21 + ^ V22D-^ (77 Vii + 5 ¥21)] F2 [V22 + e V22D-^ {rj Vu + S V22)], 
On the other hand, to have the explicit expression of the basic reproduction ratio, we need the following 
lemma whose proof is given in Appendix B. 

Lemma 3 ; Let M he a square block matrix of the following form: 

M 



Ml M2 
M3 M4 



where Mi, M2, M^ and A/4 are also square matrices. 

1. If M2 is invertihle and M2M3 - M2MiM:^^Mi = 0, then 

p{M) = max{0, p{Mi + M2M4,M^^)}. 



2. Moreover, if M2M4 = M4M2, then 



p{M) = max{0, p{Mi + M2)}. 



(15) 



(16) 



Note that A = FiiWii has the form of the matrix M defined in Lemma 3 with 
Ml = Fi [V21 + e V22D-^ (7/ Vii + 5 ¥21)] , M2 = Fi [V22 + e V22D-^ (77 V12 + 5 ¥22)] , 
Afg = F2 [V21 + i V22D-^ (7/ Vii + 5 ¥21)] and A/4 = F2 [V22 + i ¥220-^ (rj V12 + 6 V22)] ■ 
Since Fi = 13 {I — q) In and F2 = l3 q In are diagonal matrices, one has, 

A/2 A/4 = Fi [V22 + C V22D-^ (77 V12 + S V22)] F2 [V22 + C V22D-^ (77 Vi2 + S V22)] 

= Fi F2 [V22 + C V22D-^ (rj V12 + S V22)] [V22 + e V22D-^ (77 V12 + S V22)] , 
= F2 Fi [V22 + e V22D-^ (77 V12 + S V22)] [V22 + ^ V22D-^ (77 V12 + 5 ¥22)] 

= F2 [V22 + i V22D-^ (77 V12 + d V22)] Fi [V22 + ^ V22D-^ (77 Vi2 + 5 ¥22)] , 

= A/4 A/2, 

and 

A/2 A/3 - A/2 A/4 A/2" ^ Ml = A/2 A/3 - A/4 A/2 A/2~ ^ Ml , 

= A/2 A/3 -A/4 A/i, 

= ^^1 [V22 + e V22D-^ (7/ V12 + 5 ¥22)] F2 [V21 + e V22D~^ (7/ Vii + 5 ¥21)] - 

[1^22 + e "1^22^5-1 (»7 Vi2 + 5 V22)] Fi [V21 + e V22D~^ (77 Vii + (5 F21)] , 

= ^^1 F2 [V22 + e V22D-^ (77 T/12 + 5 ¥22)] [V21 + e V22D-^ (7/ Fii + (5 ¥21)] - 
F2 Fi [V22 + C ^22-D"' {v V12 + 5 ¥22)] [V21 + e 1^22^"' (?/ 1^11 + 5 V21)] , 
0. 
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(17) 



With this ill mind, since A > 0, by applying Lemma 3, Eq. (13) becomes 

= p{Fi [V2i+^V22D-^{ijVn+SV2i)] + [V22 + ^V22D-^ {r^Vr2 + SV22)]) , 

= /3p[(l-g) [V21 + ^V22D-^ {7^Vn + SV21)] +q[V22+^V22D~U7]Vr2+SV22)]] 
= /3p [(1 - q)V2i + qV22 + e V22D-^ [7? ((1 - q)Vn + qVr2) + S{{1 - q)V2i + qV22)]] , 
= /3 p [(/„ + S^V22D-^) ((1 - q) V21 +qV22)+^S V22D-' ((1 - q)Vn + qVi2)] ■ 

From the above expressions, it is evident that Vn, V12, V21, V22, i*^^ > 0. We point out that as Vu and 
V22 are irreducible and nonnegative, one has V12, V21 > 0. This imphes that A is also irreducible and non- 
ncgativc. Then, using the Perron- Frobcnius theorem [31], one can deduce that p{A) is a positive eigenvalue 
of A. Additionally, a simple calculation can prove that 

'{l-q)a{l~e)+qAE' 



[il-q)V2i+qV22] = 



a A} 



(1 - g) a (1 - 6) [aDs + b {Ae - De)] + 9 fe^B (A^ - De) 



C 



a (a - 5) Ae {Ae - De) 

= 0.5 In + h 

In+5^V22D-^ = I,,+^5iaoIn+boC){aiIn+b4C) 

= [{1+^5 0004) In+ (ao&4 + ^004 + bob4)C] 

= 06 In + be C, 

^SV22D-^ [(1 - q)Vn + qVu] = [{ae - 1) /„ + be C] [((1 - 9)^3 + gaa)/,. + ((1 - 9)^*3 + g&2)] 

= ("6 - 1) [(1 - 9)03 + qa2]In 

+ [{ae - 1) [(1 - q)b3 + qb2] + be[{l - q)a3 + qaa] + be[{l - q)b3 + 963]] C 

= 07 /„ + 67 C. 



Finally 



13 {In +S^V22D-^) [(1 - q) V21 +qV22] +^SV22D-^ [(1 - q)Vu + qVi, 

= P [(as In + b5 C){ae In + beC)+ 07 /„ + 67 C] 
= /3 [(05 ae + 07) /„ + (05 be + fos ae + &5 be + h) C] 
= asln + bs C. 
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where 

as 



{l^q)a{l~d) + qAi 



clAe 

(1 - g) a (1 - 9) [aPE + b {Ae - De)] + g b^g (^g - De) 
' " a(a-5)AB(AB-i?B) 

ae = {l + ^Saoai) 

be = S,S{aob4 + bQa4 + bobi), (18) 

a? = (ae - 1) [(1 - g)a3 + ga2], 

67 = [(ae - 1) [(1 - 9)63 + 9^2] + 66[(1- 9)03 + ga2]+66[(l- 9)63 + 9^2]], 

ag = asfle + ay, 

bs = 05 ^6 + 65 ae + 65^6 + br. 
Now, since C is a rank-one matrix that admits 1 as a unique positive eigenvalue, the greatest eigenvalue of 

the matrix is /3 [(og In + ag C] is /3 [ag + bs] and consequently, the basic reproduction ratio of system (8) is 

TZq = P [ag + 6g] , 

(19) 

= /3[(as + &5)(ae + 6e) + a7 + 67]- 

□ 

The following result is established from Theorem 2 of [30] : 

Lemma 4 ; The disease-free equilibrium Qq of system (8) is locally asymptotically stable whenever TZq < 1; 
and instable if TZq > 1 . 

Biologically speaking. Lemma 4 implies that TB can be eliminated from the community (when TZq < 1) 
if the initial sizes of the population are in the basin of attraction of the discasc-frec equilibrium Qo- 

Now, let us analyze the basic reproduction number (19). The parameter values used for numerical 
simulation are given in Table 1. 

Table 1: Description of parameters of model system 



Parameter 


Description 


Estimated value 


Source 


A 


Recruitment rate 


1001 year"^ 


[35] 


/3 


Transmission coefficient 


Variable 






Per capita naturally death rate 


0.017 year-i 


[34] 


q 


Fast route to active TB 


0.015 


[36] 


a 


Slow route to active TB 


0.0024 year-i 


[35] 





Per capital rate of effective chemoprophylaxis 


0.001 year"! 


[36] 


S 


Recovery rate of infectious 


0.7372 year-i 


[35] 


V 


Recovery rate due to chemoprophylaxis 


year~^ 


[35] 


7 


Natural recovery rate of infectious 


0.7372/4 year-i 


Assumed 


d 


Per capita disease-induced mortality rate 


0.0012 year-i 


[35] 




Relapse of recovered individuals 


0.0986 year-i 


[35] 



Figure ?? shows the effects of the transmission rate /? and the patch connectivity k on the basic reproduc- 
tion ratio 72.0 given as in Eq. (19). We have taken a metapopulation with scale- free distribution p{k) ~ 
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with (fc) = 6, k„iin = 3 and Ds = De = Dj = = 1. All other parameters are as in Table 1. The part 
above the unity of the picture corresponds to the region of the instability of the disease-free equilibrium, 
while the part below the unity of the figure represents the region for the stability of the disease-free equilib- 
rium. From this figure, one can see that TZq decreases if /3 decreases even in the case of large values of k. 
This means that if the transmission coefficient /3 is sufficiently small, TB infection could be eliminated in 
the host population even if the number of the patch connectivity k is large. However, it is difficult to control 
/?. This figure also shows that for the chosen parameter values, if the patch connectivity k does not exceed 
1.2 [k < 6), then TB can be controlled irrespective of the value of /3. The infection will equally persist for 
k > 6. 

The combined effects of the patch connectivity k and the recovery rate S on the basic reproduction 
number TZq when /? ~ 0.0017 are shown in Fig. ??. This figure suggests that the basic reproduction ratio TZq 
decreases if S increases or k decreases. Thus, the treatment of TB will have beneficial effects on infectious 
populations if the recovery rate is large. 

3.2 Analysis of the simple mass action (or density-dependent) model 

In this section, we consider the analysis of the spread of TB in metapopulation uncorrelated networks 
under the assumption of simple mass action (or density-dependent). Under these assumptions, system (8) 
can be written as 

Ps,k = A - /3 pi,kps,k - fJ-PSM - Ds (^ps,k - Jj^Ps^ ' 

/ k 

PE,k = /3(1 - q) Pi,kPs,k + /3 (1 - £,)Pi,k PR,k + 7P/,fe - [m + ^ + "(1 - (^)]PE,k - De I^Pem - 'JJ^Pe 
Pi,k = I3q Pi,kPs,k + "(1 - 0)pE,k - (/i + + 7 + S)pi,k - Di (^pi^k - '^Pi^ + ^ P^'^' 

(fc) 

(20) 

Using the same notations as in Eq. (10), system (20) can be written in the following compact form: 
f S = AI~ l3diag{I)S - {p + Ds)S + DsCS, 

E = /3{l-q)diag{I)S + /3{l-£,)dia.g{I)R + jI-[p + ri + a{l-0)+DE]E + DECE, 

(21) 

/ = l3qdiag{I)S + a{l-e)E-{p + d + -f + S + Di)I + DiCI + ^R, 

R = -f3{l-£,)dmg{I)R + jjR + SI-{p + ^ + DB)R + DRCR, 
where S, E, I, R and diag(/) are defined as in Eq. (10). 

3.2.1 Local stability of the DFE 

We give the formulae of the basic reproduction number, TZq, for the density-dependent model, using again 
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k 

PR,k ~ i)pi,k PR.k + 'n PR^k+ 5 pi^k- {P + C) PRM - Dr[ pR,k - -r^PR 



the next generation approach, developed in [30] . Then, derive bounds on TZq in term of the connectivities of 
patches. 

Using the notations in [30], the matrices F and V, for the new infections and the remaining transfers, 
are defined analogously as for the frequency-dependent model except that 

Fi = (3{1 - q)diag{S°) and F2 ^ /3qdiag{S°). 

Similarly, the techniques used in the previous subsection can be used to compute the basic reproduction 
number of the density-dependent model (20). Hence, Lemmas 1, 2 and 3 can be used to find the spectral 
radius of the following matrix: 

L = /3diag(50) ((/„ + 6iV22D-^) [(1 - q) V21 + 9^22] +^'5^22^?-^ [(1 - q)Vn + qVi2]) , 

(22) 

= 13 [asdiag{S°) + bsdiagiS")C] , 
where as and bs are defined as in Eq. (18). Thus 

TZo^^p [{as diag(5") + bs diag(5°) C] . (23) 

Since the spectral radius of L is very difficult to compute, we shall only give some properties and estimates 
of its eigenvalues owing to its specific form. 

We point out that L is a sum of a diagonal matrix ag diag(S'*') and a rang 1 matrix bs diag(5°) C. Moreover, 
the diagonal elements of og diag(S'°) arc positive and written in the increasing order. Thus, L can be 
considered as a diagonal matrix perturbed by a rank-one matrix. Now, for a general interlacing theorem of 
eigenvalues for perturbations of a diagonal matrix by rank-one matrices [33], the eigenvalues A^^ < < 
• ■ • < ^fcri ~ -^fcinax ^ iuterlacc with the eigenvalues /3 as S^^ < /? ag S^^ < ... < 13 as 5*"^ of (3 as diag(S'") 
as follows 

I3as < Xk, < l3asS'^^ < < ... < Pas S^^ < Afe„ = Afe„^^. 

Then, it follows that all the eigenvalues of L are real, simple, positive and the greatest one is Xk„-,^^ = Ti-o- 
Thus, the following inequality for 72-0 holds: 

Note that Ps.k^^^ is defined as 

Therefore, a sufficient condition for the DFE to be unstable is given by the following lemma: 
Lemma 5 // 

A(M (fc) +Ds k^ax) ^ , i^.. 

pas TTYT — —f^. — > 1, (24 

then the disease-free equilibrium of the density dependent model is unstable. 
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Condition (24) implies that TZq > 1, which is a sufficient condition for the DFE to be unstable. Rearranging 
condition (24) gives 



1 

Pas' 



(25) 



Condition (25) simply says that, if the number of individuals inhabiting those patches with highest 
connectivity in the metapopulation, for fixed values of fi, 7, De, Dj, Dfi, S, ^, q, 9, f3, d, 7, 77 and a, a 
large enough j. guarantee the instability of the disease-free equilibrium. This implies that the infection 
reaches all patches. 

Now, we prove that TZq is bounded above and below and give a sufficient condition of the instability of 
the DFE in term of the average density of patches of lowest connectivities. 

Observe that L is nonnegative (L > 0) and is an increasing function of the connectivity k. Thus 



Since S2 = 



as{mmS^k)In + bs diag(5") C 



one has 



<L</3 



as (max52) /„ + bs diag(S'") C 

k 



(3 [^L.„«8 In + bs diag(5°) C] < L < p [Sl^^^^as /„ + bs diag(5°) C] 



Then, one can deduce that 



p p K.„«8 In + 68 dmg{S°) C] < p{L) <l3p [S°,_as /„ + bg diag(5°) C] 



which implies 



/3 



' '^"min + ^8 ^ <5'" Ckk 



< P{L) < P 



•^fcmax +^8^5"° Ckk 



We have established the following lemma which give precise bounds on TZq and then yield a sufficient condition 
for the instability of the DFE in term of the average density of patches of lowest connectivities. 



Lemma 6 .■ The basic reproduction number of the density- dependent model satisfies 

.0 fcp(fc) 



/3 



kpjk) 



<TZo<P 



''sSl_+bsJ2s°k 



(26) 



The proof of this lemma is straightforward since bs diag(S'*')C is a rank one matrix, therefore, the only 

Sk ~~77T~ '^^ diagonal entries. 
. . . 

From this Lemma 6, we deduce a sufficient condition of the instability of the DFE in term of the average 
density of patches of lowest connectivities given as: 

.0 kp{k) 



Since 5*^ = Ps k' condition (27) becomes 



(k) 



> 1. 



(27) 



> 



4-^8 Y.ps,k 



kpik) 
(k) 



(28) 
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This condition (28) says that, if the number of individuals inhabiting those patches with lowest connectivity 
in the metapopulation, for fixed values of /i, 7. De, Di, Dji, (5, ^, q, 9, /3, cZ, 7, 7] and a, a large enough 
P% fcmin guarantee the instability of the disease-free equilibrium. This implies that the infection reaches all 
patches. 

In summary, it is classically known that if 7?.o < 1, then the DFE is locally stable, and if TZq > 1, then it is 
unstable. With this classic result in mind and the bounds on TZa giving by condition (24) and condition (28), 
we have established the following result giving sufficient conditions for the instability of the DFE. 

Theorem 1 .• For the model with density- dependent model (20), 
if the average density 0/ patches with highest connectivities satisfies 

> ^ (29) 

or, 

if the average density of patches with lowest connectivities satisfies 



as 



(30) 



then, the DFE is unstable. 



Model of this type demonstrates clear infection threshold. In the presence of a threshold, disease eradica- 
tion requires the reduction of the infection rate below a critical level where a stable infection-free equilibrium 
is guaranteed. In epidemiological terminology, the infection threshold may be expressed in terms of the ba- 
sic reproductive ratio 7?.o, the average number of infections produced by a single infected individual in a 
population of susceptible. From this definition, it is clear that TB infection can spread in a population only 
if 7^o > 1. 

In conclusion, crossing the threshold reduces the basic reproductive ratio TZq below unity and the infection 
is prevented from propagating. 

3.2.2 Endemic equilibrium 

Herein, we investigate the existence of an endemic equilibrium of system (21) in the special case where 
there is no re-infection after recovery (i.e. no flow from the recovered class to the latently infected class due 
to infection), but with possible relapse from the disease. 

To this end, it is more convenient to write system (21) in a more compact form. In a more compact form, 
model (21) may be written as follows: 

i: = AI - diag(B y)x +[DsC - (^t + Ds)] x, 

(31) 

il = I By){ei \ x) K-i ~Vy, 



20 



where x ^ S e M"q, y = {E, /, i?)^ e K>'o, JCi e M^" are constant vectors with 

/Ci = ( l-g,0,--- ,0 , g,0,y ,q ,o,--^- ,0/, 
IC2 = (0,l-9,0,--- ,0,0,g,0,--- ,0,0,--- ,0)^, 



/C„ - ( O,--- ,0,l-g , 0,.-.^,0,g ,0,--^- ,0/, 

is the canonical basis of M", B = [0, P /„, 0] with a n x n nuU matrix, I is defined as in Eq. (21) and V 
is the 3n x 3n constant matrix: 

'Ae Iu-DeC -7/„ 
-a (1 - 0) /„ Ai I„ -DiC 

-V -Sin Aiilr^-DRC 

We point out that the matrix —V is a Metzler matrix, that is, a matrix with all its ofi-diagonal entries 
nonncgativc [31-34]. 

With this new notations, and using the method of [30], the basic reproduction ratio (23) satisfies 



V = 



Uo= p 



(32) 



where .t" = 5" = (p^_fc)fc. 

Let Q* ~ (x*,y*) be the positive endemic equilibrium of system (31). Then, the positive endemic 
equilibrium (steady state with y > 0) can be obtained by setting the right hand side of equations in system 
(31) at zero, giving 



AI - diag(B y*) x* + [DsC-{ii + Ds)In] x* = 0, 

n 

E(e,; I By*){e, \ x*) K., -V y* = Q. 

i=l 

Multiplying the second equation of (33) by yields 

y*^Y.<e, I By*){e, \ x*)V-'lC,. 

i=l 

Using the first equation of (33), one has 

a;* = [diag(By*) -[DsC-{fi + Ds)In]]-^M. 

Then, one can deduce that 

71 

2/*=^(e, I By*){e,\ [diagiB y*) ~ [Ds C ~ {^l + Ds)In]]-'Al) V'^ IC, 
Remind that at the disease-free equilibrium, one has 



(33) 



(34) 



Al = -[DsC-{^l + Ds)In] » 0. 
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Plugging the above expression in Eq. (34) yields 

71 

y*=Y.{e, I By*){e, \ -[d\&g{B y*) ^ [Dg C ^ {ii + Ds)In]r\Ds C - {^i + Ds)In]x'')V-^ K, 

i=l 

Multiplying the above equation by B and setting z* = By* gives 

n 

z* = ^(e, I z*){e,\ -P-\z*)[Ds C - ifi + Ds)In]x°) BV-' JC,, 

where 



(35) 



P{z*) = diag(z*) -[DsC-ifi + Ds)In]- 

Wc give the explicit expression of the inverse matrix of P(z*) since we will need it later. Note that ^-'(2:*) has 

the form of the matrix R = U + XWZ given in Lemma 2 with U = diag[z* + (/.t + 1)5)1], X = [fci, ^2, . . . , fcn], 
Ds 

W = 1 and Z = ——[p(fci), 73(^2), . . . ,p(fc„)]. Then, using Lemma 2, a simple computation gives 



P^\z*) ^ diag 



1 



zl + fi + Ds 





1 


DsCdisig 




_zl+ ^ + Ds 



kp{k) 



(36) 



Now, from Eq. (35), one has 

n 

(e, I z*) =^(e, I z*)(e, I (z*)P(0) a;") ( e, | BV~^IC,), i = l,2,...,n, 
1=1 

where P(0) = — [-Ds C — (/i + Ds)/„]. From the above equation, one can deduce that 



Cj z 



I I p-i(z*)P(0)x°)^^e, I BV-'lC, 



(37) 



(38) 



Then, to find the endemic equilibrium of system (21), it suffices to find solutions of the following equation: 



where 



H{z*) 



H{z*)^l, 



^(e, I z*){e, I p-^{z*)P{0)x"){ ^ | BV-^IC, 
i=i \ j=i 

n 



(39) 



(40) 



where P~-^{z*) is defined as in Eq. (36). Note that z* are the intersection points between the curve of H(z* 
and the line 2 = 1. 

From Eq. (40), it follows that the function H{z*) satisfies 

lim H{z*) = 0, 



and 



Hm H{z*) = Y,{e., | x") ( J] e, | B V~' /C, 

^ i=i \ i=i 

We claim the following result. 
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Lemma 7 ; The inequality lim H(z*) > TZq holds. 

z'—fO 

n 

Proof: Let A = Y.{^i I B V''^ IC, ej . Then, using Eq. (32), one has 7^o = piA). Since A is 



i=l 



nonnegative matrix, if rj — ^ Aij is the sum of the j*'^ column of A, one has 



minjrj} < p{A) < max{rj}. 

j j 



If Cj denotes the canonical basis of R", I = (ei + 62 + • ■ • + e„) , using the fact that efl—1, Vi, one has 



T 

F 



n 

n 

3 I E(e. I xO)BV^-i/C,, 

i=l 
n 

E (e. I a;0)(e, | BV'^ K.,) 



With this mind, one can deduce that 



n In 

= lim 

Then, one has that 

n 

T^-o = p(^) < max{r,} < > r,-, 
which implies that lim H{z*) > TZq- This completes the proof. 

□ 

Note that we use the expression of to put emphasis on the fact that V^^ > because —V is a 
Mctzler matrix. Since lim H{z*) > TZq and lim H{z*) ~ 0, H{z*) is a positive function. Thus, positive 
solutions of Eq. (39) exist if and only if lim H{z*) > TZq > 1. From the first equation of (33), one has 

Z*— !-0 

X* = P^^{z*) AI. Since P^^{z*) is a positive definite matrix, one has x* > 0. On the other hand, since 
z* are the intersection points between the curve of H{z*) and the line z ~ 1, one has that z* > 0. Then, 
when TZq > 1, the equilibria are endemic. This means that there exists at least one endemic equilibrium of 
the model (21). Also, note that z* = By* is not a bijection (it is a onto map, but not a one to one map), 
one can conclude that the TB model with simple mass action transmission could have multiple endemic 
equilibria. However, to know the number of endemic equilibria, we need to analyze the function H{z*). We 
stress that Eq. (39) is very difficult to solve analytically due to the fact that if is a highly nonlinear function. 
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Nonetheless, one can numerically plot this curve and examine how the intersection point(s) with the line 
z = 1 change with model parameters. We have established the following theorem for the density-dependent 
model (21). 

Theorem 2 ; For the model with density- dependent model (20), if the basic reproduction number TZq > 1, 
then there exists at least one endemic equilibrium. 

3.2.3 Numerical studies 

To illustrate the various theoretical results contained in the previous section, system (20) are simulated 
using the parameter value/range in Table 1. In all simulations, the initial conditions have been chosen 
randomly. We have also taken a mctapopulation with scale-free distribution p(fc) ~ with (k) — 6 and 

kmin — 3. 

Figure ?? gives the evolution of the model (20) when (3 = 0.0001 and Ds = De = Dj ~ Dr = 1 (so that 
TZo < 1). All other parameters are as in Fig. ??. Figure ??(a) presents the prevalence curves of the model 
while, the time evolution of the number of infected individuals in each patch is depicted in Fig. ??(b). From 
these figures, it clearly appears that the disease disappears in the host population even for higher values of 
the patch connectivity. 

Figure ?? gives the evolution of the model (20) when /3 = 0.001 and Ds = De = Dj = Dji = 1 (so 
that TZf) > 1). All other parameters are as in Fig. ??. From this figure, one can observe that the disease 
persists in the host population. In addition, one can also observe that as the patch connectivity increases, 
the prevalence of the infection also increases. 

Now, let us examine the influence of the migration on the propagation of TB in the host population. 

Figure ?? presents the prevalence of the infection of model (20) in nodes of degree k of an uncorrelated 
scale-free network for different values of the migration rates. From this figure, the role of the migration rates 
Ds, De, Di and De. is remarkable. Increasing the value of the migration rates Ds, De, Dj and Dr causes 
a reduction in the prevalence of the infection. This is the only case wc have observed in which the infection 
prevalence changes non-uniformly across the mctapopulation when varying the value of a parameter. 

4 Conclusion 

In this paper, we have presented a system of differential equations of reaction-diffusion type describing 
the TB spread in heterogeneous complex metapopulations. The spatial configuration is given by the de- 
gree p{k) and the conditional probabilities P{k \ k'). For uncorrelated networks under the assumption of 
standard incidence transmission, wc have computed the disease-free equilibrium and the basic reproduction 
number. We have also shown that the disease-free equilibrium is locally asymptotically stable. Moreover, 
for uncorrelated networks and under assumption of simple mass action transmission, necessary and sufficient 
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conditions for the instability of the disease-free equihbrium for uncorrelated networks have been given in 
term of the highest and lowest connectivities of patches. We have also shown that if the basic reproduc- 
tion number TZq > 1, then the simple mass action model could have multiple endemic equilibria. Through 
numerical simulations, we found that the prevalence of the infection increases with the path connectivity. 
Also, increasing the value of the migration rates cause a reduction in the prevalence of the infection. 
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Appendix A: Proof of Lemma 1 

In this appendix, we give the proof of Lemma 1. Note that the matrix N can be written as 

'A^i N2 

'Ni 
N3 I 



N 



D 



Then, one can deduce that 
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This ends the proof. 



□ 



Appendix B: Proof of Lemma 3 

In this appendix, we give the proof of Lemma 3. To do so, we shall use the properties of the determinant. 
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(-1)" dot 



Let A the spectrum of M . Assume that M is a 2ri x 2n square matrix, then, 
det(A/-A/2„) = det 

M3 Mi-XIn_ 

AU-XIn M3 

= (-1)" det(M2) det [A/3 - (A/4 - A/„) M'^ (A/i - A/„)] , 

= (-1)" det [A/2 A/3 - A/2 A/4 A/2"^ A/i + A (A/i + A/2 A/4 M^^) - A2/„] 
If A/2 A/3 - M2 Mi A/i = 0, then 

det(A/ - A/2„) = (-A)" det [A/i + A/2 A/4 M'^ - A/„] . 

Moreover if A/2 A/4 = A/4 A/2 then 

det(A/ ~ A/2„) = (-A)" det [A/i + A/4 - A/„] . 

This ends the proof. 



□ 



Appendix C: Proof of Lemma 2 

In this Appendix, we give the proof of Lemma2. To do so, it suffices to verified that GG^^ = /„. Indeed, 
one has 

I -1 



xwzu-^x[w-^ + zu-^x] zu-\ 



In-X 
In-XW 
In - XW 



[w-^ + zu-^x] ^ + w- wzu-^x [VF-1 + zu-^x] ' 



zu-\ 



w-^ [w-^ + zu-^x] ^ -ir, + zu-^x[w-^ + zu-\x] ^ 



zu- 



[w~^ + zu~^x] [w-^ + ZU~^X] ^ - I„ 



zu~^, 



In-XW {In -In)ZU~\ 



In- 



This concludes the proof. 
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